home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 2 / BBS in a box - Trilogy II.iso / Files / Education / M / MathPad 2.15 / examples / Complex Roots next >
Encoding:
Text File  |  1993-09-19  |  5.0 KB  |  137 lines  |  [TEXT/MPAD]
  1. -- This example gives formulas for quadratic and cubic roots and uses the image command to visualize a complex function. Utility functions for complex numbers are at the end of the document.
  2.  
  3. ------------------- quadratic roots ---------------
  4. -- Algorithm for real parts of roots is from by W.H. Press, S. Teukolsky et al,"Numerical Recipes".
  5. -- Occasionally 4ac << b, so one of the roots is (erroneously) called 0.
  6. -- This formulation avoids the problem.
  7. -- implemented by David Derbes for MathPad
  8.  
  9. -- Given a quadratic of the form
  10.  
  11. --          a*x^2 + b*x + c = 0
  12.  
  13. -- with real coefficients, find the (possibly complex) roots.
  14. -- Roots are given in the form x + iy.
  15. ~
  16. sgn(x) = 1 when x >= 0,-1 otherwise
  17. D = (b*b - 4*a*c)    -- discriminant
  18. x1 = -(b + sgn(b)*sqrt(D))/(2*a) when D >= 0, -b/(2*a) otherwise
  19. x2 = c/x1 when D >= 0, -b/(2*a) otherwise
  20. y1 = 0 when D >= 0, sqrt(-D)/(2*a) otherwise
  21. y2 = -y1
  22. ~
  23.  
  24. ------------------- cubic roots -------------------
  25. -- Tartaglia's & Cardano's formulae for the roots of a cubic
  26. -- from Universal Encyclopaedia of Mathematics
  27. -- implemented by David Derbes for MathPad, 8 Sept 1993
  28. -- Given a cubic of the form
  29. --
  30. --            a0*x^3 + a1*x^2 + a2*x + a3 = 0
  31. -- 
  32. -- with real coefficients, find the (possibly complex) roots.
  33.  
  34. c1 = a1/a0  -- "normalize", i.e. make leading coefficient = 1
  35. c2 = a2/a0
  36. c3 = a3/a0
  37.  
  38. -- discriminant D; if D > 0, one real root; else three real roots
  39. -- (at most two distinct real if D = 0)
  40.  
  41. a = c2 - (c1*c1)/3.0
  42. b = (((2.0*c1*c1*c1) - (9.0*c1*c2))/27.0) + c3
  43. D = (b*b/4.0) + (a*a*a/27.0)
  44.  
  45. numRealRoots = 3 when D < 0, 1 otherwise
  46.  
  47. dHalf = sqrt(abs(D))
  48.  
  49. sgnp = -1 when ((-b/2.0) + dHalf) < 0, 1 otherwise    
  50. sgnq = -1 when ((-b/2.0) - dHalf) < 0, 1 otherwise
  51.  
  52. p = 0.0 when ((-b/2.0) + dHalf) = 0,
  53.      sgnp*(abs((-b/2.0) + dHalf))^(1.0/3.0) otherwise
  54. q = 0.0 when ((-b/2.0) - dHalf) = 0,
  55.      sgnq*(abs((-b/2.0) - dHalf))^(1.0/3.0) otherwise
  56.  
  57. s = (-b/2.0)/sqrt(-a*a*a/27.0);  theta = acos(s)
  58.  
  59. -- roots of the form x + iy
  60.  
  61. x1 = 2.0*p - (c1/3.0) when numRealRoots = 1 and abs(D) < 1.0e-10,
  62.      (p+q) - (c1/3.0) when numRealRoots = 1 and abs(D) > 1.0e-10,
  63.      2.0*sqrt(-a/3.0)*cos(theta/3.0) - (c1/3.0) otherwise
  64.  
  65. x2 = -p - (c1/3.0) when numRealRoots = 1 and abs(D) < 1.0e-10,
  66.      -(p+q)/2.0 - (c1/3.0) when numRealRoots = 1 and abs(D) > 1.0e-10,
  67.      2.0*sqrt(-a/3.0)*cos((theta/3.0) + 120) - c1/3.0 otherwise 
  68.  
  69. x3 = x2 when numRealRoots = 1,      
  70.      2.0*sqrt(-a/3.0)*cos((theta/3.0) + 240) - c1/3.0 otherwise 
  71.  
  72. y1 = 0.0  -- no matter what, must have at least one real root
  73.  
  74. y2 = (p-q)*sqrt(3.0)/2.0 when numRealRoots = 1 and abs(D) > 1.0e-10,
  75.      0.0 otherwise
  76.  
  77. y3 = -y2 when numRealRoots = 1 and abs(D) > 1.0e-10,
  78.      0.0 otherwise
  79.  
  80. root1 := {x1,y1}:;     root2 := {x2,y2}:;   root3 := {x3,y3}:
  81.  
  82. ------- Enter the coefficients for the cubic here --------------
  83. --            a0*x^3 + a1*x^2 + a2*x + a3 = 0
  84.  
  85. a0 = 35;    a1 = 15;    a2 = -5;    a3 = -20
  86.  
  87. root1:{0.76,0.00}
  88. root2:{-0.59,0.64}
  89. root3:{-0.59,-0.64}
  90.  
  91. ----------- Check the solution------------
  92. z(x) = a0*Ccube(x) + a1*Csqr(x) + a2*x + {a3,0} -- The complex cubic
  93.  
  94. -- confirm that z(x) is zero at the roots
  95. z(root1):{0.00,0.00}
  96. z(root2):{0.00,0.00}
  97. z(root3):{0.00,0.00}
  98.  
  99. -------- check the solution graphically
  100. -- define an array that samples points on the surface of:
  101. --    Z = abs(z(x)) vs. X = real(x) , Y = imaginary(x)
  102. -- This surface should dip to zero at the roots. (The sampled surface will dip near zero but in general is not sampled exactly at a root).
  103.  
  104. surf[ix,iy] =  x:=Cscale(ix,iy), Cabs(z(x)) dim[m,m]
  105.  
  106. -- The index for surf[] runs from 1 to m. Scale the index to get real and imaginary parts ranging from Xmin to Xmax and Ymin to Ymax
  107. scale(i,min,max,nsteps) = (i-.5)*(max-min)/nsteps+min
  108. Cscale(ix,iy) = {scale(ix,Xmin,Xmax,m), scale(iy,Ymin,Ymax,m)}
  109.  
  110. m=12;   -- surface is sampled on an m by m grid
  111. Xmin = -1; Xmax = 1
  112. Ymin = -1; Ymax = 1
  113. Zmin =  0; Zmax = 50
  114. image surf                  -- show image of the surface
  115. plot {root1,root2,root3}    -- show locations of roots
  116.  
  117. plot z({X,0})[1]/Zmax     -- plot the real part of z for real x
  118. -- This should be zero at real roots. On the plotted surface, real roots are located along y=0 so the real cubic plotted in this way should pass though its real roots.
  119. plot 0
  120.  
  121. --------------------------------------------------------------------
  122. ----------- utility functions for complex numbers ------
  123. -- Arrays can be used to represent complex numbers. Addition, subtraction and multiplication by a real can be done directly. The following functions implement other basic operations on complex numbers. Note that a real number r must be written as {r,0}.
  124.  
  125. Cabs(A) = sqrt(A[1]^2+A[2]^2)
  126. Cmult(A,B) = {A[1]*B[1]-A[2]*B[2],A[2]*B[1]+A[1]*B[2]}
  127. Csqr(A) = {A[1]*A[1]-A[2]*A[2],2*A[1]*A[2]}
  128. Ccube(A) = Cmult(A,Csqr(A))
  129. -- (the following were not used in this example)
  130. Cdiv(A,B)={(B[1]*A[1]+B[2]*A[2])/(B[1]^2+B[2]^2),
  131.          (B[1]*A[2]-B[2]*A[1])/(B[1]^2+B[2]^2)}
  132. Conj(A) = {A[1],-A[2]}
  133. Carg(A) = 0 when A[1]=0 and A[2]=0,
  134.          atan(A[2]/A[1]) when A[1]>=0,
  135.          atan(A[2]/A[1])+180 when A[2]>=0,
  136.          atan(A[2]/A[1])-180
  137.